#Loading librarieslibrary(sf)library(ggplot2)#Turn "s2" off to enable polygon simplificationsf_use_s2(FALSE)#Reading Spain bordersesp0 <-st_read(dsn="./data/gadm41_ESP_shp/gadm41_ESP_0.shp", quiet =TRUE)esp0<-st_simplify(esp0, dTolerance=0.01)#Note: Spain is made out of the polygons#We will focus on the largest polygons#that are traditionally association with Spain#Turning multipolygons into polygonsesp0b =st_cast(esp0,"POLYGON")#Calculating polygon areaesp0b$area<-st_area(esp0b)#Turning the units into square kmesp0b$area2<-units::set_units(esp0b$area, km^2)#Only choosing the largest polygonesp0c<-subset(esp0b, area2==max(esp0b$area2))#Making a box object around the polygonbbox <-st_bbox(esp0c)# Extracting min and max latitudes and longitudesmin_lat_y <- bbox["ymin"]max_lat_y <- bbox["ymax"]min_lon_x <- bbox["xmin"]max_lon_x <- bbox["xmax"]
#Reading wildfiresfires <-read.csv('./data/fires_2021/viirs-snpp_2021_Spain.csv')#Turning the relevant time variable into a date variable recognized by Rfires$acq_date<-as.Date(fires$acq_date)#Turning the CSV file into an SF objectfires_sf<-st_as_sf(fires, coords =c("longitude", "latitude"))#Associating the crs information for the Spain shapefile to pointsfires_sf <-st_set_crs(fires_sf, st_crs(esp0)) # 4326 is the EPSG code for WGS84#Only selecting fires of type 0#See the variable definitions:#0 = presumed vegetation fire#1 = active volcano#2 = other static land source#3 = offshore detection (includes all detections over water)fires_sf<-subset(fires_sf, type==0)# Extract month informationfires_sf$month <-format(fires_sf$acq_date, "%m")#Reading only Augustfires_aug<-subset(fires_sf, acq_date>"2021-07-31"& acq_date<"2021-09-01")
#Defining an error for better map extenterror<-0.05fires_aug$bright_ti4ggplot() +geom_sf(data = esp0, color ="black", fill="black") +geom_sf(data = fires_aug, aes(color = bright_ti4, size = bright_ti4), alpha =0.7) +scale_color_gradient(low ="yellow", high ="darkred", name ="bright_ti4") +scale_size_continuous(range =c(0.5, 2), name ="bright_ti4") +#We still need to define the extent: fires also happen outside the large polygon#associated with Spaincoord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error)) +ggtitle("Map with Dots Colored and Sized by bright_ti4") +guides(color =guide_legend(title.position ="top"),size =guide_legend(title.position ="top"))
Making the map
We can finally make the map:
Entire Year
We can also produce the following map
Entire Year
This is done using facet_wrap
ggplot() +geom_sf(data = esp0, color ="black", fill="black") +geom_sf(data = fires_sf, aes(color = bright_ti4, size = bright_ti4), alpha =0.7) +scale_color_gradient(low ="yellow", high ="darkred", name ="bright_ti4") +scale_size_continuous(range =c(0.5, 2), name ="bright_ti4") +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error)) +ggtitle("Map with Dots Colored and Sized by bright_ti4") +guides(color =guide_legend(title.position ="top"),size =guide_legend(title.position ="top"))+facet_wrap(~month, ncol =4) # Adjust the number of columns as needed
Entire Year
This is done using facet_wrap
ggplot() +geom_sf(data = esp0, color ="black", fill="black") +geom_sf(data = fires_sf, aes(color = bright_ti4, size = bright_ti4), alpha =0.7) +scale_color_gradient(low ="yellow", high ="darkred", name ="bright_ti4") +scale_size_continuous(range =c(0.5, 2), name ="bright_ti4") +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error)) +ggtitle("Map with Dots Colored and Sized by bright_ti4") +guides(color =guide_legend(title.position ="top"),size =guide_legend(title.position ="top"))+facet_wrap(~month, ncol =4) # Adjust the number of columns as needed
Entire Year
Or we can produce an animated map:
Entire Year
Or we can produce an animated map:
library("gganimate")#Extracting coordinatesfires_sf$lat_x<-st_coordinates(fires_sf)[,"X"]fires_sf$lon_y<-st_coordinates(fires_sf)[,"Y"]#Dropping the geometryfires<-st_drop_geometry(fires_sf)#Mappingfig<-ggplot() +geom_sf(data = esp0, color ="black", fill="black")+#We need to use geom_point NOT geom_sf because of library compatibility problemsgeom_point(data = fires, aes(x= lat_x, y = lon_y, color = bright_ti4, size = bright_ti4), alpha =0.9)+#Note that animated geom_sf with points does not work because of some library incompatibility#geom_sf(data = fires_sf, aes(color = bright_ti4, size = bright_ti4), alpha = 0.7) +scale_color_gradient(low ="yellow", high ="darkred", name ="bright_ti4") +scale_size_continuous(range =c(0.1, 4), name ="bright_ti4") +guides(color =guide_legend(title ="bright_ti4"),size =guide_legend(title ="bright_ti4"))+coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error))+#animate argumentstransition_time(acq_date) +ggtitle('Date: {frame_time}')#defining animation options: speed of transition and reolutionanimate(fig, fps=2, res =60)
Spatial Join
spatial join is one of the most common operations in spatial analysis.
In a spatial join we are “attaching” attributes from one layer to another, based on their spatial relations.
For example, we will match the fires to Spanish provinces.
Spatial Join
Spatial Join
We can now try to merge the points to level-2 subdivisions: gadm41_ESP_2
Spatial Join
Place the folder with the shape files in the relevant week folder
#Reading wildfiresfires <-read.csv('./data/fires_2021/viirs-snpp_2021_Spain.csv')#Turning the acq_date variable into a date variablefires$acq_date<-as.Date(fires$acq_date)#Make sf object from pointsfires<-st_as_sf(fires, coords =c("longitude", "latitude"))fires <-st_set_crs(fires, st_crs(esp0)) # 4326 is the EPSG code for WGS84# Extract month and year informationfires$month <-format(fires$acq_date, "%m")fires$year <-format(fires$acq_date, "%Y")#Simplifying our datafires<-subset(fires, type==0)#Reading only Augustfires_aug<-subset(fires, acq_date>"2021-07-31"& acq_date<"2021-09-01")
Spatial Join
We will try to count the number of fires by province
Spatial Join - Step 1
And this is how we produce this graph:
Step 1:
library(dplyr)library(gridExtra)#Selecting only one variableesp2b<-subset(esp2, select =c("NAME_2"))#Selecting only bright_ti5. #This represents I-5 Channel brightness temperature of the fire pixel measured in Kelvin.fires_augb<-subset(fires_aug, select =c("bright_ti5"))#Perfoming the spatial joinxjoin<-st_join(fires_augb, esp2b)
The unit is fire (i.e. the invidual event). The unit is NOT the district.
This is why Avila appears over 1,000 times
Thus, we need to count how many fires we have per district.
Spatial Join - Step 1
#Dropping the geometry from the object (i.e. turning it into a regular dataframe)#We need to drop the geometry to aggregate by groupxjoin<-st_drop_geometry(xjoin)#Summing by polygonssum_per_province <- xjoin %>%group_by(NAME_2) %>%summarize(sum_fires =sum(fire_no, na.rm =TRUE))#Performing a left join to the shapefilesesp1c<-left_join(esp2b, sum_per_province, by =c("NAME_2"="NAME_2"))#Replacing NAs with 0esp1c$sum_fires[is.na(esp1c$sum_fires)]<-0#Calculating logs to eliminate extreme valuesesp1c$log_fires<-log(esp1c$sum_fires+1)
Spatial Join - Step 2
And this is how we produce this graph:
Step 2:
#Calculating SD and Meanfires_sd <-sd(esp1c$log_fires, na.rm=T)fires_mean <-mean(esp1c$log_fires, na.rm=T)#Creating different stepsstep1<-min(esp1c$log_fires, na.rm=T)step2<-fires_mean - fires_sdstep3<-fires_mean -0.5*fires_sdstep4<-fires_meanstep5<-fires_mean +0.5*fires_sdstep6<-fires_mean + fires_sdstep7<-max(esp1c$log_fires, na.rm=T)#Putting all the steps in a listclint <-c(step1, step2, step3, step4, step5, step6, step7)clint
#Reading Spain bordersesp0 <-st_read(dsn="./data/gadm41_ESP_shp/gadm41_ESP_0.shp", quiet =TRUE)esp0<-st_simplify(esp0, dTolerance=2000)#Simplifying country polygonesp0b =st_cast(esp0,"POLYGON")#Calculating polygon areaesp0b$area<-st_area(esp0b)#Turning the units into square kmesp0b$area2<-units::set_units(esp0b$area, km^2)#Only choosing the largest polygonesp0c<-subset(esp0b, area2==max(esp0b$area2))#Making a box object around the polygonbbox <-st_bbox(esp0c)# Extracting min and max latitudes and longitudesmin_lat_y <- bbox["ymin"]max_lat_y <- bbox["ymax"]min_lon_x <- bbox["xmin"]max_lon_x <- bbox["xmax"]
```{r , fig.width = 6.5, fig.height=2.4}library(gridExtra)#Arranging the Map and Histogram together#in grid.arrange there are 4 slots in a row. you need to indicate which slot will be taken#by fig1 vs. fig2. In our case, fig1 takes the first 3 slots#fig2 takes the 4th slot.grid.arrange(p1,p2, layout_matrix=rbind(c(1,1,1,2)))```
Spatial Join - Step 6
And this is how we produce this graph:
Step 6:
Spatial Join - Step 6
To save the map that you made still look good, you should use the fillowing dimensions:
st_buffer calculates buffers - polygons encompassing all points up to the specified distance from a given geometry.
We can for example calculate 250km buffers around Avila and Valencia.
To be precise about our distances, we will use EPSG codes
# Transform your data to a suitable projected CRS that uses meters as units# For example, you can use a UTM zone that covers your area of interest# Here, I'm assuming your data is around Avila, Spain, so I'm using UTM zone 30N (EPSG:2062)# See https://epsg.io/ and look for Spain# You should choose an appropriate UTM zone for your data. Note that otherwise, this would be in degreesavila_ctr_sp <-st_transform(avila_ctr, crs =2062)valencia_ctr_sp <-st_transform(valencia_ctr, crs =2062)#Creating the buffer at 250km. Note that unit is in meters: st_crs(avila_ctr_sp) - LENGTHUNIT under Datumavila_100 =st_buffer(avila_ctr_sp, dist =250000)valencia_100 =st_buffer(valencia_ctr_sp, dist =250000)#Going back to the original crs in degreesavila_100 <-st_transform(avila_100, crs =st_crs(esp1c))valencia_100 <-st_transform(valencia_100, crs =st_crs(esp1c))
Buffers
ggplot()+geom_sf(data=esp1c)+geom_sf(data=avila_ctr, color ="red", size=3)+geom_sf(data=valencia_ctr, color ="red", size=3)+geom_sf(data=avila_100, color ="red", fill=NA)+geom_sf(data=valencia_100, color ="red", fill=NA)+#geom_sf(data=l, color = "red")+theme_void() +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error))
Buffers
Another type of operation is buffer
st_buffer calculates buffers - polygons encompassing all points up to the specified distance from a given geometry.
We can for example calculate 250km buffers around Avila and Valencia:
Geometry Difference from Pairs
There are four geometry-generating functions that operate on pairs of input geometries:
st_intersection
st_difference
st_sym_difference
st_union
Geometry Difference from Pairs
Geometry Difference from Pairs
Geometry Difference from Pairs
Geometry Difference from Pairs
Geometry Difference from Pairs
Geometry Difference from Pairs: union
We might be interested in the total area that is within 250km from both Anvila and Valencia
We can do this in the following way:
union_area<-st_union(avila_100, valencia_100)ggplot()+geom_sf(data=esp1c)+geom_sf(data=avila_ctr, color ="red", size=3)+geom_sf(data=valencia_ctr, color ="red", size=3)+geom_sf(data=union_area, color ="red", fill=NA)+theme_void() +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error))
Geometry Difference from Pairs: union
We might be interested in the total area that is within 250km from both Anvila and Valencia
We can do this in the following way:
Geometry Difference from Pairs: union
Remember that st_union can also be applied to an individual layer,
This returns a dissolved union of all geometries:
Vector layer aggregation
Sometimes, we may not want to dissolve all features into a single geometry
Sometimes, we may want to dissolve by an attribute
Vector layer aggregation
Let us assume that we did not have the polygons for provinces in Spain.
A common operation is to calculate density of lines within polygons
For example, we would be interested in calculating density of railways within provinces
#Step1: Reading the polylinesrails <-st_read(dsn="./data/europe-railways-shape/railways.shp", quiet =TRUE)#Step2: Calculating the length of rails in mlengths <-st_length(rails)#Step3: Uniting the linesmerged_line <-st_union(rails)#Step4: Summing up the lengths of the original linestotal_length <-sum(lengths)#Step5: Creating a new sf object with the merged line and total lengthresult_sf <-st_sf(geometry = merged_line, total_length = total_length)ggplot()+geom_sf(data=esp1c)+geom_sf(data=result_sf, color ="red")+theme_void() +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error))
Vector layer aggregation
A common operation is to calculate density of lines within polygons
For example, we would be interested in calculating density of railways within provinces
Vector layer aggregation
#Step1: Performing the intersection between rails and polygonsintersection <-st_intersection(rails, esp1c)#Step2: Calculating the length of the intersecting polylinesintersection$length <-st_length(intersection)#Step3: Saving only relevant varsintersection2<-subset(intersection, select =c("NAME_2", "length"))#Step4: Dropping the geometryintersection2<-st_drop_geometry(intersection2)#Step5: Dividing the length by 1000 to obtain kmintersection2$length<-as.numeric(intersection2$length)/1000#Step6: Performing the left join (left = polygons; right = lines)esp1c<-left_join(esp2b, intersection2, by =c("NAME_2"="NAME_2"))#Mappingggplot(esp1c, aes(fill=length)) +geom_sf() +geom_sf_text(aes(label = NAME_2), colour ="white", size=2)+theme_void() +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error))+scale_fill_viridis_c(option ="viridis")
Vector layer aggregation
Conclusion
We have covered a veriety of sf functions:
join data to a vector by location
agregate points by polygon
perform geometric calculations: areas, change units